Optimisation of planing¶

This page was generated from notebook Models/Basic_France_model/Planing_optimisation/case_planing_step_by_step_learning.ipynb with nbconvert

If you're lost, you can go back to

  • one node model README it gives planing tools and consumption models
  • two node models README for step by step learning with two nodes (France and Germany).
  • Mode complete 7 node model README for a more complete European model to do prospective analysis

Table of Contents¶

  • 1. Introduction
  • 2. First basic problem
    • 2.1. Math and first step with pyomo for solving the problem
    • 2.2. Analysing results : lagrange multipliers
  • 3. Extensions of this planing problem
    • 3.1. Linear temporal coupling with ramp constraints
  • 4. Case with storage
  • 5. System with 100% of renewable
  • 6. 2035, 4 million EV, 30 TWh H2, 40 GW nuke DSM

1. Introduction ¶

This document contains a description of optimisation tools for electric system planing simulation with only one node. Having one node is a strong approximation but it is usefull to learn how the tool works.

If you don't known anything about these optimisatino problem, you might prefer to start with operation optimisation in file optim-Operation.ipynb, you can have a look at the whole notebook's results through the web page here. In files BasicFunctionalities/input-XXX.ipynb you can learn to understand input data (consumption, availability).

(Text below is almost the same as for operation)

This document will gives a chance to understand

  • how to do the optimisation of planing of an electric system under the hourly operation constraint
  • understand the mathematical formulation of the optimisation problem
  • learn to analyse de results of the optimisation and in particular the Lagrange multipliers
  • get in touch with pyomo (a python package to write optimisation problems)

It proposes to enter the subject by increasing progressively the number of variables and constraints in the optimisation problem, hence moving toward more realism through the document, introducing:

  • ramp constraints that implies a simple temporal coupling
  • spatially indexed variables and congestion constraints that implies a simple spatial coupling.
  • storage constraints that implies a temporal coupling

It relies on different test cases that allow to

  • consider different production means (nuclear, CCG, solar, onshore wind power, offshore wind power, hydro, curtailement of consumption, storage)
  • consider different meteorological years for France
  • consider different countries in the multi-zone case (France, Germany, GB, Spain)

If, after reading this file, you want to build your own pyomo model you can start by having a look at GetElectricSystemModel_PlaningSingleNode function in f_planingModels.py

In [1]:
#region importation of modules
import os
import sys
if os.path.basename(os.getcwd())=='Planing_optimisation':
    sys.path.append('../../../')
import sys
if sys.platform != 'win32':
    myhost = os.uname()[1]
else : myhost = ""
if (myhost=="jupyter-sop"):
    ## for https://jupyter-sop.mines-paristech.fr/ users, you need to
    #  (1) run the following in a terminal
    if (os.system("/opt/mosek/9.2/tools/platform/linux64x86/bin/lmgrd -c /opt/mosek/9.2/tools/platform/linux64x86/bin/mosek.lic -l lmgrd.log")==0):
        os.system("/opt/mosek/9.2/tools/platform/linux64x86/bin/lmutil lmstat -c 27007@127.0.0.1 -a")
    #  (2) definition of license
    os.environ["MOSEKLM_LICENSE_FILE"] = '@jupyter-sop'

import numpy as np
import pandas as pd
import csv
#import docplex
import datetime
import copy
import plotly.graph_objects as go
import matplotlib.pyplot as plt
from sklearn import linear_model
import sys


from functions.f_graphicalTools import *
from functions.f_consumptionModels import *
from Models.Basic_France_models.Planing_optimisation.f_planingModels import *
# Change this if you have other solvers obtained here
## https://ampl.com/products/solvers/open-source/
## for eduction this site provides also several professional solvers, that are more efficient than e.g. cbc
#endregion
In [2]:
#region Solver and data location definition
InputConsumptionFolder='../Consumption/Data/'
InputProductionFolder='../Production/Data/'
InputPlaningFolder='Data/'
GraphicalResultsFolder="GraphicalResults/"
InputFolder='Data/input/'

solver= 'mosek' ## no need for solverpath with mosek.
BaseSolverPath='/Users/robin.girard/Documents/Code/Packages/solvers/ampl_macosx64'
sys.path.append(BaseSolverPath)

solvers= ['gurobi','knitro','cbc'] # 'glpk' is too slow 'cplex' and 'xpress' do not work
solverpath= {}
for solver in solvers : solverpath[solver]=BaseSolverPath+'/'+solver
cplexPATH='/Applications/CPLEX_Studio1210/cplex/bin/x86-64_osx'
sys.path.append(cplexPATH)
solverpath['cplex']=cplexPATH+"/"+"cplex"
solver = 'mosek'
#endregion

2. First basic problem ¶

2.1. Math and first step with pyomo for solving problem¶

Before you start with the math, you should

  • be expert in optimisation of operation, your can learn things in optim-Operation.ipynb
  • be familiar with annualized cost of electricity (variable and fixed cost) and levelised cost of energy. You can read my post here on lcoe (in french, if there are people that would be willing to contribute to a translation that would be great !!) and also this one on the role of power and enery in system cost.
\begin{align} &\text{Cost function }& &\min_{x,\bar{x}} \sum_i \left ( \beta_i\bar{x_i}+ \sum_t\pi_i x_{it}\right ) \;\;\; & & \pi_i \text{ marginal cost, }\beta_i \text{fixed annualized cost}\\ &\text{Power limit } & &\text{ s.t.} \;\; 0 \leq x_{it}\leq a_{it} \bar{x_i} & &\bar{x_i} \text{ installed power, } a_{it} \text{ availability}\\ &\text{Meet demand } & & \sum_i x_{it} \geq C_t && C_t \text{ Consumption}\\ \end{align}
In [3]:
#region I - Simple single area : loading parameters
year=2013
#### reading areaConsumption availabilityFactor and TechParameters CSV files
areaConsumption = pd.read_csv(InputConsumptionFolder+'areaConsumption'+str(year)+'_FR.csv',sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date"])
availabilityFactor = pd.read_csv(InputProductionFolder+'availabilityFactor'+str(year)+'_FR.csv',sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date","TECHNOLOGIES"])
TechParameters = pd.read_csv(InputPlaningFolder+'Planing-Simple_TECHNOLOGIES.csv',sep=',',decimal='.',skiprows=0,comment="#").set_index(["TECHNOLOGIES"])
TechParameters.head()
#### Selection of subset
Selected_TECHNOLOGIES=['OldNuke','CCG'] #you can add technologies here
availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
TechParameters=TechParameters.loc[Selected_TECHNOLOGIES,:]
#TechParameters.loc[TechParameters.TECHNOLOGIES=="OldNuke",'maxCapacity']=63000 ## Limit to actual installed capacity
#endregion

/var/folders/zt/7cz1y9md79l_h08bbqymt4w9z8xlw7/T/ipykernel_78518/463447673.py:10: FutureWarning: The behavior of indexing on a MultiIndex with a nested sequence of labels is deprecated and will change in a future version. `series.loc[label, sequence]` will raise if any members of 'sequence' or not present in the index's second level. To retain the old behavior, use `series.index.isin(sequence, level=1)`
  availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
In [4]:
#region I - Simple single area  : Solving and loading results
model = GetElectricSystemModel_PlaningSingleNode(Parameters={"areaConsumption"      :   areaConsumption,
                                                   "availabilityFactor"   :   availabilityFactor,
                                                   "TechParameters"       :   TechParameters})

if solver in solverpath :  opt = SolverFactory(solver,executable=solverpath[solver])
else : opt = SolverFactory(solver)
results=opt.solve(model)
## result analysis
Variables=getVariables_panda_indexed(model)
print(extractCosts(Variables))
print(extractEnergyCapacity(Variables))

#pour avoir la production en KWh de chaque moyen de prod chaque heure
production_df=Variables['energy'].pivot(index="Date",columns='TECHNOLOGIES', values='energy')
### Check sum Prod = Consumption
Delta=(production_df.sum(axis=1) - areaConsumption.areaConsumption);
abs(Delta).max()
#endregion

              Capacity_Milliards_euros  Energy_Milliards_euros
TECHNOLOGIES                                                  
CCG                           1.262412               -0.000000
OldNuke                      13.894902                3.444102
              Capacity_GW  Production_TWh
TECHNOLOGIES                             
CCG                  10.8         0.00000
OldNuke              56.7       492.01462
Out[4]:
0.0

2.2 Analysing results : lagrange multipliers ¶

Verify that the sum of market prices allows all actors to cover fixed and marginal cost. do they earn more ? why ?

In [5]:
#region I - Simple single area  : visualisation and lagrange multipliers
### representation des résultats
fig=MyStackedPlotly(y_df=production_df,Conso = areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename=GraphicalResultsFolder+'file.html') ## offline
fig.show()

#### lagrange multipliers
Constraints= getConstraintsDual_panda(model)

# Analyse energyCtr
energyCtrDual=Constraints['energyCtr']; energyCtrDual['energyCtr']=energyCtrDual['energyCtr']
energyCtrDual
round(energyCtrDual.energyCtr,2).unique()

# Analyse CapacityCtr
#CapacityCtrDual=Constraints['CapacityCtr'].pivot(index="Date",columns='TECHNOLOGIES', values='CapacityCtr')*1000000;
#round(CapacityCtrDual,2)
#round(CapacityCtrDual.OldNuke,2).unique() ## if you increase by Delta the installed capacity of nuke you decrease by xxx the cost when nuke is not sufficient
#round(CapacityCtrDual.CCG,2).unique() ## increasing the capacity of CCG as no effect on prices
#endregion
Out[5]:
array([7.])

3. Extensions of this planing problem ¶

3.1. Linear temporal coupling with ramp constraints ¶

In the this section, we will increase the complexity of the problem given in Section 2 and add : dependency on area z (country),a congestion constraint, ramp constraints.

\begin{align} &\text{Cost function }& &\min_{x} \sum_z \left ( \beta_{iz}\bar{x_{iz}} + \sum_t \sum_i \pi_{iz} x_{itz} \right ) \;\;\; & & \pi_{iz} \text{ marginal cost}\\ &\text{Power limit } & &\text{ s.t.} \;\; 0 \leq x_{itz}\leq a_{itz} \bar{x_{iz}} & &\bar{x_{iz}} \text{ installed power, } a_{itz} \text{ availability}\\ &\text{Meet demand } & & \sum_i x_{itz} \geq C_{tz} && C_{tz} \text{ Consumption}\\ &\text{Stock limit } & &\sum_t x_{it}\leq E_i && E_i=\bar{x_i}*N_i \text{ Energy capacity limit}\\ &\text{ramp limit } & &rc^-_i *x_{it}\leq x_{it}-x_{i(t+1)}\leq rc^+_i *x_{it} && rc^+_i rc^-_i\text{ ramp limit}\\ \end{align}

In [6]:
#region II - Ramp Single area : loading parameters loading parameterscase with ramp constraints
year=2013
Selected_TECHNOLOGIES=['OldNuke', 'CCG',"curtailment"] #you'll add 'Solar' after
#### reading CSV files
areaConsumption = pd.read_csv(InputConsumptionFolder+'areaConsumption'+str(year)+'_FR.csv',
                                sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date"])
availabilityFactor = pd.read_csv(InputProductionFolder+'availabilityFactor'+str(year)+'_FR.csv',
                                sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date","TECHNOLOGIES"])
TechParameters = pd.read_csv(InputPlaningFolder+'Planing-RAMP1BIS_TECHNOLOGIES.csv',sep=',',decimal='.',skiprows=0,comment="#").set_index(["TECHNOLOGIES"])

#### Selection of subset
availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
TechParameters=TechParameters.loc[Selected_TECHNOLOGIES,:]
TechParameters.loc["OldNuke",'RampConstraintMoins']=0.01 ## a bit strong to put in light the effect
TechParameters.loc["OldNuke",'RampConstraintPlus']=0.02 ## a bit strong to put in light the effect
#endregion

/var/folders/zt/7cz1y9md79l_h08bbqymt4w9z8xlw7/T/ipykernel_78518/253438452.py:12: FutureWarning:

The behavior of indexing on a MultiIndex with a nested sequence of labels is deprecated and will change in a future version. `series.loc[label, sequence]` will raise if any members of 'sequence' or not present in the index's second level. To retain the old behavior, use `series.index.isin(sequence, level=1)`
In [7]:
#region II - Ramp Single area : solving and loading results
model = GetElectricSystemModel_PlaningSingleNode(Parameters={"areaConsumption"      :   areaConsumption,
                                                   "availabilityFactor"   :   availabilityFactor,
                                                   "TechParameters"       :   TechParameters})
opt = SolverFactory(solver)
results=opt.solve(model)
Variables=getVariables_panda_indexed(model)
print(extractCosts(Variables))
print(extractEnergyCapacity(Variables))

#pour avoir la production en KWh de chaque moyen de prod chaque heure
production_df=Variables['energy'].pivot(index="Date",columns='TECHNOLOGIES', values='energy')
### Check sum Prod = Consumption
Delta=(production_df.sum(axis=1) - areaConsumption.areaConsumption);
abs(Delta).max()
#endregion

              Capacity_Milliards_euros  Energy_Milliards_euros
TECHNOLOGIES                                                  
curtailment                    0.00000              687.829674
CCG                            1.40268                2.709080
OldNuke                        9.80240                3.007976
              Capacity_GW  Production_TWh
TECHNOLOGIES                             
curtailment           9.0       22.927656
CCG                  12.0       39.376163
OldNuke              40.0      429.710802
Out[7]:
1.4551915228366852e-11
In [8]:
#region II - Ramp Single area : visualisation and lagrange multipliers
fig=MyStackedPlotly(y_df=production_df,Conso = areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename=GraphicalResultsFolder+'file.html') ## offline
fig.show()

#### lagrange multipliers
Constraints= getConstraintsDual_panda(model)

# Analyse energyCtr
energyCtrDual=Constraints['energyCtr']; energyCtrDual['energyCtr']=energyCtrDual['energyCtr']*1000000
#energyCtrDual
#round(energyCtrDual.energyCtr,2).unique()

# Analyse CapacityCtr
#CapacityCtrDual=Constraints['maxCapacityCtr'].pivot(index="Date",columns='TECHNOLOGIES', values='CapacityCtr')*1000000;
#round(CapacityCtrDual,2)
#round(CapacityCtrDual.OldNuke,2).unique() ## if you increase by Delta the installed capacity of nuke you decrease by xxx the cost when nuke is not sufficient
#round(CapacityCtrDual.CCG,2).unique() ## increasing the capacity of CCG as no effect on prices
#endregion

4. Case with storage ¶

4.1. Optimisation of a storage market participation ¶

This only change is that this time, storage parameters $P_{max}$ and $C_{max}$ are decision variables.

In [9]:
#region III Ramp+Storage single area : loading parameters
Zones="FR"
year=2013

Selected_TECHNOLOGIES=['OldNuke','WindOnShore', 'CCG',"curtailment",'HydroRiver', 'HydroReservoir',"Solar"] ## try adding 'HydroRiver', 'HydroReservoir'

#### reading CSV files
areaConsumption = pd.read_csv(InputConsumptionFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',
                                sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date"])
availabilityFactor = pd.read_csv(InputProductionFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',
                                sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date","TECHNOLOGIES"])
TechParameters = pd.read_csv(InputPlaningFolder+'Planing-RAMP1BIS_TECHNOLOGIES.csv',
                             sep=',',decimal='.',skiprows=0,comment="#").set_index(["TECHNOLOGIES"])
StorageParameters = pd.read_csv(InputPlaningFolder+'Planing-RAMP1_STOCK_TECHNO.csv',sep=',',decimal='.',skiprows=0).set_index(["STOCK_TECHNO"])

#### Selection of subset
availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
TechParameters=TechParameters.loc[Selected_TECHNOLOGIES,:]
#TechParameters.loc["CCG",'capacity']=100000 ## margin to make everything work
TechParameters.loc["CCG",'maxCapacity']=70000
TechParameters.loc["OldNuke",'maxCapacity']=35000
TechParameters.loc["OldNuke",'RampConstraintMoins']=0.03 ## a bit strong to put in light the effect
TechParameters.loc["OldNuke",'RampConstraintPlus']=0.03 ## a bit strong to put in light the effect
#endregion

/var/folders/zt/7cz1y9md79l_h08bbqymt4w9z8xlw7/T/ipykernel_78518/2552567369.py:17: FutureWarning:

The behavior of indexing on a MultiIndex with a nested sequence of labels is deprecated and will change in a future version. `series.loc[label, sequence]` will raise if any members of 'sequence' or not present in the index's second level. To retain the old behavior, use `series.index.isin(sequence, level=1)`
In [10]:
#region III Ramp+Storage single area : solving and loading results
model = GetElectricSystemModel_PlaningSingleNode_withStorage(Parameters={"areaConsumption"      :   areaConsumption,
                                                   "availabilityFactor"   :   availabilityFactor,
                                                   "TechParameters"       :   TechParameters,
                                                 "StorageParameters": StorageParameters})

if solver in solverpath :  opt = SolverFactory(solver,executable=solverpath[solver])
else : opt = SolverFactory(solver)
results=opt.solve(model)
Variables = getVariables_panda_indexed(model)
Constraints = getConstraintsDual_panda(model)

production_df=Variables['energy'].pivot(index="Date",columns='TECHNOLOGIES', values='energy')
production_df.loc[:,'Storage'] = Variables['storageOut'].pivot(index='Date',columns='STOCK_TECHNO',values='storageOut').sum(axis=1)-Variables['storageIn'].pivot(index='Date',columns='STOCK_TECHNO',values='storageIn').sum(axis=1) ### put storage in the production time series
production_df.sum(axis=0)/10**6 ### energies produites TWh
production_df[production_df>0].sum(axis=0)/10**6 ### energies produites TWh
production_df.max(axis=0)/1000 ### Pmax en GW

fig=MyStackedPlotly(y_df=production_df, Conso=areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
#endregion

WARNING: Implicitly replacing the Component attribute stockLevel_ini
    (type=<class 'pyomo.core.base.param.IndexedParam'>) on block unknown with
    a new Component (type=<class 'pyomo.core.base.var.IndexedVar'>). This is
    usually indicative of a modelling error. To avoid this warning, use
    block.del_component() and block.add_component().

5. System with 100% of renewable ¶

This case is a one node case, not realistic to analyse the future of the electric system, but usefull to learn to use the tool. If you want to use a more complete case, you can use the 7 node EU model

This simple system has no interconnection, no flexible demand. Flexibility is mainly done with storage and gaz.

In [11]:
#region IV Case Storage + CCG + PV + Wind + hydro (Ramp+Storage single area) : loading parameters
Zones="FR"
year=2013

Selected_TECHNOLOGIES=['CCG', 'WindOnShore','WindOffShore','Solar',"curtailment",'HydroRiver', 'HydroReservoir']

#### reading CSV files
areaConsumption = pd.read_csv(InputConsumptionFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',
                                sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date"])
availabilityFactor = pd.read_csv(InputProductionFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',
                                sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date","TECHNOLOGIES"])
TechParameters = pd.read_csv(InputPlaningFolder+'Planing-RAMP1BIS_TECHNOLOGIES.csv',
                             sep=',',decimal='.',skiprows=0,comment="#").set_index(["TECHNOLOGIES"])
StorageParameters = pd.read_csv(InputPlaningFolder+'Planing-RAMP1_STOCK_TECHNO.csv',
                                sep=',',decimal='.',skiprows=0).set_index(["STOCK_TECHNO"])
#### Selection of subset
availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
TechParameters=TechParameters.loc[Selected_TECHNOLOGIES,:]
#TechParameters.loc["CCG",'capacity']=100000 ## margin to make everything work
TechParameters.loc["CCG",'energyCost']=100
TechParameters.loc["CCG",'maxCapacity']=50000
TechParameters.loc["WindOnShore",'capacityCost']=120000 #€/MW/year - investment+O&M fixed cost
TechParameters.loc["Solar",'capacityCost']=65000 #€/MW/year
TechParameters.loc["CCG",'RampConstraintMoins']=0.4 ## a bit strong to put in light the effect
TechParameters.loc["CCG",'RampConstraintPlus']=0.4 ## a bit strong to put in light the effect
StorageParameters.loc["Battery1","p_max"]=10000 # this is not optimized - batteries
StorageParameters.loc["Battery2","p_max"]=7000 # this is not optimized - Pumped HS
StorageParameters.loc["Battery2","c_max"]=StorageParameters.loc["Battery2","p_max"]*20 # this is not optimized 20h of Pumped HS
#endregion

/var/folders/zt/7cz1y9md79l_h08bbqymt4w9z8xlw7/T/ipykernel_78518/2233090117.py:17: FutureWarning:

The behavior of indexing on a MultiIndex with a nested sequence of labels is deprecated and will change in a future version. `series.loc[label, sequence]` will raise if any members of 'sequence' or not present in the index's second level. To retain the old behavior, use `series.index.isin(sequence, level=1)`
In [12]:
#region IV Case Storage + CCG + PV + Wind + hydro  (Ramp+Storage single area) : solving and loading results
model = GetElectricSystemModel_PlaningSingleNode_withStorage(Parameters={"areaConsumption"      :   areaConsumption,
                                                   "availabilityFactor"   :   availabilityFactor,
                                                   "TechParameters"       :   TechParameters,
                                                   "StorageParameters"   : StorageParameters})

if solver in solverpath :  opt = SolverFactory(solver,executable=solverpath[solver])
else : opt = SolverFactory(solver)
results=opt.solve(model)
Variables = getVariables_panda_indexed(model)
Constraints = getConstraintsDual_panda(model)

production_df=Variables['energy'].pivot(index="Date",columns='TECHNOLOGIES', values='energy')
production_df.loc[:,'Storage'] = Variables['storageOut'].pivot(index='Date',columns='STOCK_TECHNO',values='storageOut').sum(axis=1)-Variables['storageIn'].pivot(index='Date',columns='STOCK_TECHNO',values='storageIn').sum(axis=1) ### put storage in the production time series
production_df.sum(axis=0)/10**6 ### energies produites TWh
production_df[production_df>0].sum(axis=0)/10**6 ### energies produites TWh
production_df.max(axis=0)/1000 ### Pmax en GW

fig=MyStackedPlotly(y_df=production_df, Conso=areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
#endregion

WARNING: Implicitly replacing the Component attribute stockLevel_ini
    (type=<class 'pyomo.core.base.param.IndexedParam'>) on block unknown with
    a new Component (type=<class 'pyomo.core.base.var.IndexedVar'>). This is
    usually indicative of a modelling error. To avoid this warning, use
    block.del_component() and block.add_component().

6. 2035, 4 million EV, 30 TWh H2, 40 GW nuke DSM ¶

This case is a one node case, not realistic to analyse the future of the electric system, but usefull to learn to use the tool. If you want to use a more complete case, you can use the 7 node EU model

In [13]:
#region V - Simple single area +4 million EV +  demande side management +30TWh H2: loading parameters
Zones="FR" ; year=2013
#### reading areaConsumption availabilityFactor and TechParameters CSV files
#areaConsumption = pd.read_csv(InputConsumptionFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date"])

TemperatureThreshold = 15
ConsoTempe_df=pd.read_csv(InputConsumptionFolder+'ConsumptionTemperature_1996TO2019_FR.csv',parse_dates=['Date']).\
    set_index(["Date"])[str(year)]
ConsoTempe_df=ConsoTempe_df[~ConsoTempe_df.index.duplicated(keep='first')]
(ConsoTempeYear_decomposed_df,Thermosensibilite)=Decomposeconso(ConsoTempe_df,TemperatureThreshold=TemperatureThreshold)


#obtaining industry-metal consumption
#  & x["type"] == "Ind" & x["UsageDetail"] == "Process").\
Profile_df_sans_chauffage=pd.read_csv(InputConsumptionFolder+"ConsumptionDetailedProfiles.csv").\
    rename(columns={'heures':'Heure',"WeekDay":"Jour"}).\
    replace({"Jour" :{"Sat": "Samedi" , "Week":"Semaine"  , "Sun": "Dimanche"}}). \
    query('UsagesGroupe != "Chauffage"'). \
    assign(is_steel=lambda x: x["Nature"].isin(["MineraiMetal"])).\
    set_index(["Mois", "Heure",'Nature', 'type',"is_steel", 'UsagesGroupe', 'UsageDetail', "Jour"]).\
    groupby(["Mois","Jour","Heure","type","is_steel"]).sum().\
    merge(add_day_month_hour(df=ConsoTempeYear_decomposed_df,semaine_simplifie=True,French=True,to_index=True),
          how="outer",left_index=True,right_index=True).reset_index().set_index("Date")[["type","is_steel","Conso"]]. \
    pivot_table(index="Date", columns=["type","is_steel"], values='Conso')
Profile_df_sans_chauffage.columns = ["Autre","Ind_sans_acier","Ind_acier","Residentiel","Tertiaire"]

Profile_df_sans_chauffage=Profile_df_sans_chauffage.loc[:,Profile_df_sans_chauffage.sum(axis=0)>0]
Profile_df_n=Profile_df_sans_chauffage.div(Profile_df_sans_chauffage.sum(axis=1), axis=0) ### normalisation par 1 et multiplication
for col in Profile_df_sans_chauffage.columns:
    Profile_df_sans_chauffage[col]=Profile_df_n[col]*ConsoTempeYear_decomposed_df["NTS_C"]

steel_consumption=Profile_df_sans_chauffage.loc[:,"Ind_acier"]
steel_consumption.max()
steel_consumption[steel_consumption.isna()]=110
steel_consumption.isna().sum()
# if you want to change thermal sensitivity + add electric vehicle

VEProfile_df=pd.read_csv(InputConsumptionFolder+'EVModel.csv', sep=';')
NbVE=10 # millions
ev_consumption = NbVE*Profile2Consumption(Profile_df=VEProfile_df,Temperature_df = ConsoTempe_df.loc[str(year)][['Temperature']])[['Consumption']]

h2_Energy = 30000## H2 volume in GWh/year
h2_Energy_flat_consumption = ev_consumption.Consumption*0+h2_Energy/8760
to_flexible_consumption=pd.DataFrame({'to_flex_consumption': steel_consumption,'FLEX_CONSUM' : 'Steel'}).reset_index().set_index(['Date','FLEX_CONSUM']).\
    append(pd.DataFrame({'to_flex_consumption': ev_consumption.Consumption,'FLEX_CONSUM' : 'EV'}).reset_index().set_index(['Date','FLEX_CONSUM'])).\
    append(pd.DataFrame({'to_flex_consumption': h2_Energy_flat_consumption,'FLEX_CONSUM' : 'H2'}).reset_index().set_index(['Date','FLEX_CONSUM']))

availabilityFactor = pd.read_csv(InputProductionFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date","TECHNOLOGIES"])



TechParameters = pd.read_csv(InputPlaningFolder+'Planing-RAMP1BIS_TECHNOLOGIES.csv',sep=',',decimal='.',skiprows=0,comment="#").set_index(["TECHNOLOGIES"])
StorageParameters = pd.read_csv(InputPlaningFolder + 'Planing-RAMP1_STOCK_TECHNO.csv', sep=',', decimal='.',
                                skiprows=0).set_index(["STOCK_TECHNO"])
ConsoParameters = pd.read_csv(InputPlaningFolder + "Planing-Conso-FLEX_CONSUM.csv", sep=";").set_index(["FLEX_CONSUM"])
ConsoParameters_ = ConsoParameters.join(
    to_flexible_consumption.groupby("FLEX_CONSUM").max().rename(columns={"to_flexible_consumption": "max_power"}))

Selected_TECHNOLOGIES=['OldNuke','CCG','TAC', 'WindOnShore', 'WindOffShore','HydroReservoir','HydroRiver','Solar','curtailment']#you can add technologies here
availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
TechParameters=TechParameters.loc[Selected_TECHNOLOGIES,:]

TechParameters.loc["CCG",'energyCost']=100
TechParameters.loc["CCG",'maxCapacity']=50000
TechParameters.loc["WindOnShore",'capacityCost']=120000 #€/MW/year - investment+O&M fixed cost
TechParameters.loc["Solar",'capacityCost']=65000 #€/MW/year
TechParameters.loc["CCG",'RampConstraintMoins']=0.4 ## a bit strong to put in light the effect
TechParameters.loc["CCG",'RampConstraintPlus']=0.4 ## a bit strong to put in light the effect
StorageParameters.loc["Battery1","p_max"]=10000 # this is not optimized - batteries
StorageParameters.loc["Battery2","p_max"]=7000 # this is not optimized - Pumped HS
StorageParameters.loc["Battery2","c_max"]=StorageParameters.loc["Battery2","p_max"]*20 # this is not optimized 20h of Pumped HS


areaConsumption=pd.DataFrame(ConsoTempeYear_decomposed_df.loc[:,"Consumption"]-steel_consumption,columns=["areaConsumption"])




def labour_ratio_cost(df):  # higher labour costs at night
    if df.hour in range(7, 17):
        return 1
    elif df.hour in range(17, 23):
        return 1.5
    else:
        return 2


labour_ratio = pd.DataFrame()
labour_ratio["Date"] = areaConsumption.index.get_level_values('Date')
labour_ratio["FLEX_CONSUM"] = "Steel"
labour_ratio["labour_ratio"] = labour_ratio["Date"].apply(labour_ratio_cost)
labour_ratio.set_index(["Date","FLEX_CONSUM"], inplace=True)
#model.labour_ratio = Param(model.Date, initialize=labour_ratio.squeeze().to_dict())


# endregion

/var/folders/zt/7cz1y9md79l_h08bbqymt4w9z8xlw7/T/ipykernel_78518/812183875.py:7: FutureWarning:

Indexing a DataFrame with a datetimelike index using a single string to slice the rows, like `frame[string]`, is deprecated and will be removed in a future version. Use `frame.loc[string]` instead.

/var/folders/zt/7cz1y9md79l_h08bbqymt4w9z8xlw7/T/ipykernel_78518/812183875.py:45: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

/var/folders/zt/7cz1y9md79l_h08bbqymt4w9z8xlw7/T/ipykernel_78518/812183875.py:46: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

/var/folders/zt/7cz1y9md79l_h08bbqymt4w9z8xlw7/T/ipykernel_78518/812183875.py:60: FutureWarning:

The behavior of indexing on a MultiIndex with a nested sequence of labels is deprecated and will change in a future version. `series.loc[label, sequence]` will raise if any members of 'sequence' or not present in the index's second level. To retain the old behavior, use `series.index.isin(sequence, level=1)`
In [14]:
#region V - Simple single area +4 million EV +  demande side management +30TWh H2: solving and loading results

# €/kW/an coût fixe additionnel pour un GW d'électrolyseur en plus en supposant que l'on construit
model =  Model_SingleNode_online_flex(Parameters={"areaConsumption"      :   areaConsumption,
                                           "availabilityFactor"   :   availabilityFactor,
                                           "TechParameters"       :   TechParameters,
                                           "StorageParameters"   : StorageParameters,
                                           "to_flexible_consumption" : to_flexible_consumption,
                                           "ConsoParameters" : ConsoParameters_,
                                           "labour_ratio": labour_ratio})


if solver in solverpath :  opt = SolverFactory(solver,executable=solverpath[solver])
else : opt = SolverFactory(solver)
results=opt.solve(model)
Variables = getVariables_panda_indexed(model)
Variables.keys()
Variables['increased_max_power'] ## on a ajouté X GW à ce qui existait.
print(extractCosts(Variables))
print(extractEnergyCapacity(Variables))
Constraints = getConstraintsDual_panda(model)

production_df=Variables['energy'].pivot(index="Date",columns='TECHNOLOGIES', values='energy')
production_df.loc[:,'Storage'] = Variables['storageOut'].pivot(index='Date',columns='STOCK_TECHNO',values='storageOut').sum(axis=1)-Variables['storageIn'].pivot(index='Date',columns='STOCK_TECHNO',values='storageIn').sum(axis=1) ### put storage in the production time series
production_df.sum(axis=0)/10**6 ### energies produites TWh
production_df[production_df>0].sum(axis=0)/10**6 ### energies produites TWh
production_df.max(axis=0)/1000 ### Pmax en GW


fig=MyStackedPlotly(y_df=production_df, Conso=areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
#endregion

/Users/robin.girard/Documents/Code/Etude/Energy-Alternatives-Planing/Models/Basic_France_models/Planing_optimisation/../../../Models/Basic_France_models/Planing_optimisation/f_planingModels.py:362: FutureWarning:

pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.


                Capacity_Milliards_euros  Energy_Milliards_euros
TECHNOLOGIES                                                    
HydroRiver                      0.807500                0.000000
Solar                           2.246639                0.000000
OldNuke                         9.802400                1.772419
WindOnShore                     1.609200                0.000000
HydroReservoir                  0.565250                0.000000
curtailment                     0.000000               -0.000000
CCG                             0.000000               -0.000000
TAC                             2.556545                4.499245
WindOffShore                    4.728668                0.000000
                Capacity_GW  Production_TWh
TECHNOLOGIES                               
HydroRiver        10.000000       48.773636
Solar             34.563681       40.784792
OldNuke           40.000000      253.202748
WindOnShore       13.410000       27.350294
HydroReservoir     7.000000       14.700000
curtailment     1000.000000        0.000000
CCG                0.000000        0.000000
TAC               25.514422       53.562441
WindOffShore      22.423501       80.520562